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Abstract. We study numerically domain growth and interface fluctuations in one- 
and two-dimensional lattice systems composed of four species that interact in a cyclic 
way. Particle mobility is implemented through exchanges of particles located on 
neighboring lattice sites. For the chain we find that the details of the domain growth 
strongly depend on the mobility, with a higher mobility yielding a larger effective 
domain growth exponent. In two space dimensions, when also exchanges between 
mutually neutral particles are possible, both domain growth and interface fluctuations 
display universal regimes that are independent of the predation and exchange rates. 
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1. Introduction 

Since the seminal works of Lotka pQ and Volterra [2], many in-depth studies of model 
systems yielded important insights into food webs and other ecological systems. Using 
the methods of non-linear and statistical physics, these investigations allowed to make 
notable progress in our understanding of biodiversity, species coexistence, and species 
extinction [31 HI El EJ [7] . Recent years saw a flurry of studies of systems composed of 
multiple species that interact in a cyclic way. Most of these studies focused on the case 



of three species 0El[ini[IIl[I2[ia[IH[Ii2[ig[I^ 

Ell EHl ES EQl EH E21 [331 EH ESI EH1 EH EHl [39], a special situation where every species 




interact with every other species. Only rather few papers, however, dealt with more 
realistic cases where a given species interacts with only a subgroup of all species living 



in the same ecological environment (6J El El EHJ SQj SH S21 SHJ S31 SHI S3, SHI EDl EH 
[521 [531 EH [55J [56] . 




The case of four species interacting in a cyclic way is characterized by the 
formation of two alliances of neutral species, i.e. species that are not in a predator- 
prey relationship. In the well-mixed situation [BTl l52l 155] a variety of orbits in 
configuration space are obtained in mean field approximation. Stochastic evolution 
of finite systems shows a good agreement with the mean field results, but also 
reveals interesting additional aspects, most notably those related to various extinction 
scenarios. Restricting themselves to predation events, Frachebourg et al. [HI [9] studied 
one-dimensional coarsening and the related algebraic domain growth. Segregation 
and the formation of defensive alliances in two dimensions have been studied in 
[3H IH1 113 SHI S?l SH] for a variety of interaction schemes and/or different realizations 
of particle mobility. 

The studies of cyclically competing species have allowed to address a variety of 
problems, ranging from biodiversity and species extinction to the formation of space- 
time pattern. In the present work we address generic properties of pattern formation 
during a coarsening process that involves multiple species. For this we study domain 
growth and interface fluctuations in one- and two-dimensional lattice systems when four 
species compete cyclically. 

Calling the four species A, B, C, and D, we prepare fully occupied lattices of N sites 
by placing for each of the four species particles on N/4 randomly selected lattice sites. 
Here A is the predator of B, whereas B preys on C, C on D, and D on A. We update 
the system by randomly selecting a pair of neighboring sites. If the individuals on the 
selected two sites are in a predator-prey relationship, then one of the following three 
events can occur: predation takes place where the prey is replaced by a predator, the 
two individuals swap places, or nothing happens. This can be written in the following 
way: 





A + B A + A 
A + B + A 
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and so on. Here, k^B and sab are the predation and swapping rates between an A and 
a B particle, with Icab + $ab < 1- If, on the other hand, the two selected individuals are 
neutral, as it is the case for A and C particles, for example, then they can exchange places 
with rate s n [17]. The results for domain growth and interface fluctuations discussed in 
the following have been obtained for cases where all the predation rates are identical, 
k = kAB = ksc — k>cD = k D A- In addition, we also set all the swapping rates between 
predator-prey pairs to be equal, s = sab = $bc — s cd = s^a- In this way, we avoid 
having a bias in favor of one of the alliances. 

The remainder of the paper is organized in the following way. In the next Section 
we discuss the impact of mobility on the coarsening process in one dimension. When 
exchanges between neutral particles are forbidden, an effective algebraic growth law is 
encountered where the exponent increases with increasing swapping rate. If we also 
allow for exchanges between neutral pairs, the domains grow exponentially fast. The 
two-dimensional case is the subject of Section 3. For that case we also find that adding 
swappings between neutral pairs changes qualitatively the domain growth. Indeed, in 
case of identical predation and exchange rates for all predator-prey pairs, the system ends 
up in a steady state where all species coexist. However, when we allow for exchanges 
between neutral pairs, coarsening sets in where the domains only contain mutually 
neutral partners. In variance with the one-dimensional case the exponent of the algebraic 
growth in two dimensions is found to be independent of the values k and s of the 
predation and swapping rates. Section 4 focuses on the interface fluctuations between 
domains composed of partner-pairs. We thereby find that the interface fluctuations 
belong to the Edwards- Wilkinson universality class [57]. Finally, in Section 5 we discuss 
our findings and conclude. 

2. Coarsening in one dimension 

In low dimensions stochastic effects have a remarkable impact on the properties of non- 
equilibrium systems, as evidenced for example by the non-trivial nature of many non- 
equilibrium phase transitions taking place in one dimension [58J. Coarsening processes 
in one dimension are usually substantially different from those taking place at higher 
dimensions. In a series of papers Frachebourg et al. [H HI H2] studied the coarsening 
process in one- dimensional systems composed of multiple species that interact in a cyclic 
way. Restricting themselves to immobile particles, they found algebraically growing 
domains both for three and four species, but with different values of the growth exponent: 
3/4 for three species and 1/3 for four species. For larger numbers of species, the system 
rapidly settles into a blocked configuration where neighboring pairs are non-interacting. 

Recent studies of three species systems have shown the importance of particle 
mobility when cyclic dominance prevails [T5| [T8| [T6| [21] . Thus in two dimensions 
mobility induces qualitative changes in cases where the total number of particles is not 
conserved. Whereas for low mobility biodiversity prevails, medium values of the mobility 
lead to species extinction. Further increasing the mobility finally yields an enhanced 
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mixing paired with the (re)emergence of species coexistence. In one dimension one finds 
for the case of conserved particle number that coarsening takes place at low mobility, 
paired with a power-law increase of the average domain size. For high swapping rates, 
however, complex space-time patterns emerge that ultimately yield a non-equilibrium 
steady state where all three species coexist. 

In order to get a first impression on how mobility changes one-dimensional 
coarsening for four species we can have a look at the space-time diagrams shown in 
Figure [TJ These diagrams have been obtained for small systems composed of N — 2500 
sites, with periodic boundary conditions. As usual, time is measured in Monte Carlo 
Steps (MCS), with one MCS corresponding to N proposed updates. In all three cases 
shown in that figure, we have set k + s — 1. The left diagram shows the typical time 
evolution of the system for small swapping rates. As for the case without swapping 
[9] we see rapidly the emergence of partner-pair domains (green and yellow versus red 
and blue) that are very stable against attacks of the competing alliance. Consider for 
example the case where a yellow and a blue domain are in contact. As the B particles 
(yellow) are the preys of the A particles (blue), the blue domain expands into the yellow 
domain until this front meets some D particles (green domains). Now, the roles are 
changed, and the blue domain decreases as the A particles are replaced by D particles. 
This continues until the expanding green domain encounters a red domain, and so on. 
This rather complicated coarsening process keeps going on until only one alliance fills 
the complete system. Increasing the swapping rate leaves this scenario unchanged (see 
middle diagram), but coarsening seems to proceed in an accelerated way. Indeed, a 
higher swapping rate allows the prey to survive with a finite probability the passage 
through small domains of its predators, so that the protection due to the pairing of 
neutral partners becomes less effective. For very high swapping rate, see the right 
diagram, it takes some time before first domains form. Once these domains are formed, 
the coarsening proceeds very rapidly. 

This discussion can be made more quantitative by looking at the average domain 
size, (A), see Figure [2} for various cases with k + s — 1. We thereby define as average 
domain size the average length of the segments composed by only one species. For 
small swapping rates the domain growth proceeds as for the immobile case, yielding an 
algebraic relationship (A) ~ t s , where S = 1/3 is the exactly known value when s = [9], 
see the dashed line. Increasing the value of s has two notable effects. First, the increase 
at early times gets slower and slower, thus revealing that the enhanced mobility makes 
it more difficult for the initial domains to form. Once ordered segments have formed, 
however, the growth proceeds much faster than for the s = case. In order to see that 
we show in the inset the time dependence of the effective exponent 

5 eff (t) = [ln((A)(10t)) - ln«A>(t))] / [ln(10t) - \n(t)] . (1) 

We see that this time-dependent effective exponent decreases for increasing times. In 
fact, the existing data leave open the possibility that asymptotically the value 1/3 is 
recovered for extremely long times. 
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Figure 1. Space-time diagrams (time increases from top to bottom) for one- 
dimensional systems with three different predation and exchange rates: k = 0.8, 
s = 0.2, s n = (left), k = 0.1, s = 0.9, s„ = (middle), k = 0.01, s = 0.99, 
Sn = (right). The system size is N = 2500, the total simulation time is t max = 10000 
MCS. Periodic boundary conditions are used. Blue: A particles, yellow: B particles, 
red: C particles, green: D particles. 

This behavior, which is due to the formation of neutral alliances, is very different 
to what is observed in the three species case, see Figures 1 and 2 in [21]. Indeed, as in 
the three species case every species interact with every other, one does not have regions 
formed by domains of mutually neutral species. As a result increasing the swapping 
rate in that case does not yield an acceleration of the coarsening, and the domain 
growth exponent remains unchanged. It is only for very large swapping rates that a 
notable change sets in, yielding a non-equilibrium steady state with constant average 
domain size, as the system is then well mixed [21]. We expect the scenario observed in 
the present paper for the four species case to be very generic for systems with cyclic 
dominance composed of multiple species, as long as different species can form neutral 
alliances to fight off their predators. 

Following [17] we can also allow for exchanges between neutral pairs, i.e. swappings 
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Figure 2. Average domain size as a function of time for different swapping rates s, 
with k + s = 1. The systems contain 36000 lattice sites. The dashed line indicates the 
exponent 1/3 known exactly for the case of immobile particles. The inset shows the 
time-dependence of the effective exponent ([1]). The shown data have been obtained 
after averaging over typically 3000 different realizations of the noise. 

of neighboring A and C particles or B and D particles with rate s n . In that case we 
observe, see Figure [3j that the domain growth proceeds in an explosive way, yielding 
a typical domain size that increases exponentially with time. Obviously, exchanges 
between neutral partners weakens the lines of defense that result from alternating 
domains formed by the partners, see Figure [H which drastically increases the probability 
of dissolving smaller domains. 

3. Coarsening in two dimensions 

A variety of studies have clarified some of the properties encountered when putting 
M > 3 species with cyclic dominance on a square lattice. Early on Frachebourg and 
Krapivsky [12] found that in the absence of particle mobility systems with more than 14 
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Figure 3. Allowing for swapping between neutral pairs yields an exponential increase 
of the domain size, as shown in this linear-log plot. Systems with 36000 sites are 
considered. 

species end up in a frozen state where all particles are surrounded by neutral partners. 
For the case we are interested in, namely that of four species, the system settles into 
a steady state where all four species coexist. Allowing for mobile particles through the 
diffusion via unoccupied sites yields a qualitative change [33]. For small densities of 
empty sites the behavior is similar to that of immobile particles, and all four species 
coexist. For higher densities, however, a coarsening regime is revealed where partner- 
pairs are competing against each other, yielding the formation of domains that contain 
only neutral species. A similar transition between a regime of species coexistence and a 
regime where species extinction proceeds through a coarsening process is observed for 
a completely filled lattice when exchanges between neutral partners take place [37] . No 
exchanges between predators and preys were allowed in that work. 

In the following we study the coarsening process taking place in our system when 
exchanges between neutral partners are allowed in addition to interactions (predation 
and exchange) between predators and preys. We thereby focus on two different sets of 
rates, namely (1) k = 0.8 and s = s n = 0.2 and (2) k = 0.2 and s = s n = 0.8, similar 
results being obtained for other values of predation and swapping rates. We stress 
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that, in agreement with [27], coarsening does not take place in absence of exchanges 
between neutral partners, but instead the system then ends up in a stationary state 
with coexitence. 







Figure 4. Configurations during coarsening in two dimensions. Top row: k = 0.8 and 
s = s n = 0.2, bottom row: k — 0.2 and s = s n = 0.8. The left column shows the 
systems 100 MCS after initialization, the snapshots in the right column are taken after 
500 MCS. The system size is 400 x 400, and periodic boundary conditions are used. 



Let us again start by having a look at some configurations, see Figure HI where 
the two different cases are displayed. The first row shows two snapshots, one after 100 
MCS and one after 500 MCS, for k = 0.8 and s = s n = 0.2. As expected, one observes 
the formation of domains that are composed of mutually neutral species, namely A and 
C (blue and red) or B and D (white and green). These domains then coarsen until, 
eventually, only one alliance will prevail in our finite system. For the case k = 0.2 and 
s = s n = 0.8, shown in the second row, domains are not yet well formed after 100 
MCS. At a later time we again have coarsening, where smaller domains give place to 
larger domains, but the boundaries of the domains are very fuzzy. These effects are 
of course due to the enhanced swapping rate and are similar to those observed in the 
one-dimensional case, see Figure [U 

Due to the fact that the boundaries of the domains are not very sharp, we determine 
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Figure 5. Average correlation length as a function of time for the same two sets of 
predation and exchange rates as those used in Figure HI The dashed line indicates the 
slope 1/2. The system size is 800 x 800, and the data result from an average over 8000 
independent runs. 



the typical length in our system through a study of the correlation function, 

C(t, r) = [{m(t, r H(o> o)) - (m(t, r)> (m(0, 0)>] . (2) 

i 

Here Ui(t, r) is an occupation number that is 1 when at time t site r is occupied by a 
particle of species i and otherwise. (• ■ ■) indicates an average over initial states and over 
the realizations of the noise. We then define the correlation length X(t) as the distance 
at which C (t, r) drops to one third of the value it has at r = 0. The result of this analysis 
is shown in Figure El At short times the correlation length is much smaller for the larger 
swapping rate, see the black line in the figure, in accordance with our observation that 
in that case it takes longer for domains to form. Once these domains are formed, the 
coarsening undergoes a regime of rapid growth before entering the asymptotic regime. 
This behavior is in fact very similar to what happens in the d = 1 case, see Figure 
EJ However, whereas in one dimension the asymptotic regime is characterized by an 
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algebraic growth with mobility dependent exponents, in two dimensions the asymptotic 
growth is the same for the studied cases, with (A(i)) ~ t 1//2 . This universal exponent 
1/2, which has also been observed in other variants of the four species model on the 
square lattice [6], is expected from the Allen-Cahn law for curvature-driven coarsening 
[591 El EI]. 



4. Interface fluctuations 

Even so the same growth law prevails in two dimensions for different reaction rates, the 
boundaries between the different domains can look very differently, see the configurations 
shown in Figure HI Whereas for small swapping rates these boundaries are very sharp, 
they get increasingly fuzzier when increasing this exchange rate. It is therefore a priori 
not clear that interfaces for different rates have the same properties. 

In order to elucidate this we prepare our system such that in one half we have only A 
and C particles, whereas the other half is formed by the competing alliance composed of 
B and D particles. We start with a sharp, straight interface, using reflective boundaries 
in the direction perpendicular to the interface and periodic boundary conditions parallel 
to the interface. We then update the system using the same rules as before and monitor 
the roughening of the interface, see Figure [6] for two examples. 

Interface fluctuations are usually studied in a quantitative way through the interface 
width. For large swapping rates, however, the interface is very diffuse, due to the 
many exchanges taking place at the boundaries between the two alliances. In order to 
determine the local position of the interface, we assign the value +1 to every A and C 
particle and the value —1 to every B and D particle and then proceed as in [62] (see 
also [63J). Considering rectangular systems of L x H spins, we determine for each row 
j the value I that minimizes the sum 

v(t) = j:iS l , J -p(z-£)} 2 . (3) 

i=l 

Here, Sij = ±1, depending on which alliance occupies site (i, j), and p{u) is a step 
function, with p = 1 for u < and p = — 1 for u > [62J. The local position of the 
interface at row j, £(j), is then given by the value of i that minimizes ([3]). With this 
coarse-grained interface profile at hand, we compute the time-dependent interface width 



W(t) 



1 JL, s _ X 2 



\ i =1 



H 

where £ = -n J2 £(J) is the mean position of the fluctuating interface at time t. 

i=i 

Figure [7] shows the time evolution of the interface width for the case k = 0.8, 
s = 0.2, and s n = 0.2. We observe three different regimes, as expected for the width of an 
interface that initially was given by a straight line: uncorrelated fluctuations prevail at 
early times, followed by a correlated fluctuation regime, before the fluctuations saturate 
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Figure 6. Snapshots of two systems where each half of the system is initially occupied 
by one alliance. Left: k — 0.5, s = 0.5, and s n — 0.5; right: k = 0.2, s = 0.8, and 
s n = 0.8. The system size is 200 x 400, the configurations are obtained after 980 MCS. 

at a level that depends on the length H of the interface. The two last regimes can be 
summarized by the Family- Vicsek scaling relation [6U [65] 

W(t,H) = H a f{t/H a ^) , (5) 

with the growth exponent /3 and the roughness exponent a. It follows from this relation 
that the width should grow as v in the correlated regime and that the value of saturation 
should scale as H a with the interface length. Analyzing our data, we find in the two 
studied cases that (3 = 1/4 (this is indicated by the dashed line in Figure [7]) and a = 1/2, 
in agreement with the values expected for the Edwards- Wilkinson universality class [57J . 

5. Conclusion 

When four species compete cyclically, the formation of two alliances of mutually neutral 
partners is observed. As a result, species extinction, when it takes place, proceeds 
through a coarsening process where domains formed by the different alliances expand 
until one alliance fills the whole system. 

In this work we studied the properties of this coarsening process, both in one and 
two dimensions. It is known from earlier studies [SJ [9] that in the case of immobile 




Figure 7. Time evolution of the interface width for the case k = 0.8, s = 0.2, and 
s„ = 0.2. Rectangular systems with L x H sites are considered, where L — 200. The 
dashed line indicates the expected power law for correlated fluctuations belonging to 
the Edwards- Wilkinson universality class. The data result from averaging over at least 
2000 independent realizations of the noise. 



particles the four species form domains in one dimension that grow algebraically in time 
with an exponent 1/3. Mobile particles render the formation of initial domains more 
difficult. Once these domains are formed, however, the growth proceeds very fast. If only 
predator-prey pairs are allowed to exchange sites, then we obtain an effective power law 
exponent whose value increases with the swapping rate. If in addition neutral pairs can 
swap places, then the domains grow exponentially fast. In two dimensions the system 
settles in a steady state characterized by the coexistence of all four species if no exchanges 
between neutral partners take place. Allowing for these exchanges, though, yields the 
formation of domains formed by partner-pairs, followed by a coarsening process. Even 
so the boundaries between the domains can change their character, being sharp for 
small exchange rates, but fuzzy for larger rates, the domain growth exponent takes on 
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the value 1/2, as expected for curvature-driven coarsening. In addition, after coarse- 
graining the interface, we find that the interface fluctuations are characterized by the 
exponents of the Edwards- Wikinson universality class. 

Obviously, the properties of four species that interact in a cyclic way are very 
different from those encountered for the case of three species. Three species form in fact 
a very special situation, as here every species interact with every other species. In four 
and more species, however, there are mutually neutral species. As a consequence one 
observes for an even number of species the formation of alliances of non-interacting 
species, thereby enhancing the efficiency of the different species to fight off their 
predators. As this is a very general mechanism, we expect our results for four species to 
be representative also for other systems, as long as the species can arrange themselves 
in two alliances composed exclusively by partners that are mutually neutral. 
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